  /***************************************************************************
  *   Copyright (C) 2007 by Alexandre Masselot,,,   *
  *   alex@alex-laptop   *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
  *   it under the terms of the GNU Library General Public License as       *
  *   published by the Free Software Foundation; either version 2 of the    *
  *   License, or (at your option) any later version.                       *
  *                                                                         *
  *   This program is distributed in the hope that it will be useful,       *
  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  *   GNU General Public License for more details.                          *
  *                                                                         *
  *   You should have received a copy of the GNU Library General Public     *
  *   License along with this program; if not, write to the                 *
  *   Free Software Foundation, Inc.,                                       *
  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
  ***************************************************************************/
#ifndef SPECTRUMALIGNMENT_H
#define SPECTRUMALIGNMENT_H
#include "math.h"
#include <iostream>

#include "spectrum.h"
#include "deltamass.h"

  namespace bol{
  /**
   * Align 2 spectra a build an index to see to which index of spectrumB matches element i of spectrumA (-1 if no peak is matched)
   *
   * For a given peak in A, we align with the highest intensity (ot the closest if intensites are equals)
   *
   * We of course assume spectra are sorted

  @author Alexandre Masselot,,, <alex@alex-laptop>
   */
  template<class A, class B> class SpectraAlignment{
    const Spectrum<A*>* m_spectrumA;
    const Spectrum<B*>* m_spectrumB;
    vector<int> m_index;

    public:
        /**
     * builds an alignement (which will anyway be read only after this stage)
     * @param spA the source spectra
     * @param spB the spectra to be matched against
     * @param tol the tolerance to allow for a peak/peak match
         */
      SpectraAlignment<A, B>(const Spectrum<A*>* const spA, const Spectrum<B*>* const spB, const DeltaMass* tol):m_spectrumA(spA), m_spectrumB(spB){
        align(tol);
      }
      SpectraAlignment<A, B>():m_spectrumA(NULL), m_spectrumB(NULL){
      }
/*
        ~SpectraAlignment<A, B>(){
        cerr<<"remove spextraalignement"<<endl;
      }
*/
      
      void align(const DeltaMass* tol){
        int nA=m_spectrumA->size();
        int nB=m_spectrumB->size();
        m_index.assign(nA, -1);
        //m_index.clear();
        int iA=0;
        double massA=(nA>0)?(*m_spectrumA)[iA]->moz():1.e99;
        int iB=0;
        double massB=(nB>0)?(*m_spectrumB)[iB]->moz():1.e99;
        double minfB, msupB;
        tol->range(massB, minfB, msupB);
        while(1){
          while(massA<minfB){
            m_index[iA]=-1;
            //m_index.push_back(-1);
            iA++;
            if(iA>=nA)
              break;
            massA=(*m_spectrumA)[iA]->moz();
          }
          if(iA>=nA)
            break;

          while(massA>msupB){
            iB++;
            if(iB>=nB)
              break;
            massB=(*m_spectrumB)[iB]->moz();
            tol->range(massB, minfB, msupB);
          }
          if(iB>=nB)
            break;
        
          if(tol->isInRange(massA, massB)){
             //locate the most intense (or the closest) B peak in the neighborhood)
            int imax=iB;
            double intmax=(*m_spectrumB)[iB]->intensity();
            double deltamin=fabs(massA-massB);
            int i=iB;
            double minfA, msupA;
            tol->range(massA, minfA, msupA);
            for(int i=iB; i<nB && (*m_spectrumB)[i]->moz()<=msupA; i++){
              if((*m_spectrumB)[i]->intensity()>intmax || ((*m_spectrumB)[i]->intensity()==intmax && fabs(massA-(*m_spectrumB)[i]->moz())<deltamin)){
                imax=i;
                intmax=(*m_spectrumB)[i]->intensity();
                deltamin=fabs(massA-(*m_spectrumB)[i]->moz());
              }
            }
            iB=imax;
            m_index[iA]=iB;
            //m_index.push_back(iB);
            iA++;
            if(iA>=nA)
              break;
            massA=(*m_spectrumA)[iA]->moz();
            massB=(*m_spectrumB)[iB]->moz();
            tol->range(massB, minfB, msupB);
          }
        }
      }

      vector< int > index() const {
        return m_index;
      }

      const Spectrum<A*>* const spectrumA() const {
        return m_spectrumA;
      }

      void spectrumA(const Spectrum<A*>* const spA){
        m_spectrumA=spA;
      }
      void spectrumB(const Spectrum<B*>* const spB){
        m_spectrumB=spB;
      }

      const Spectrum<B*>* const spectrumB() const {
        return m_spectrumB;
      }
      template<class U, class V> friend ostream& operator<<(ostream& out, const SpectraAlignment<U, V>& spalign);

  };
  template<class U, class V> ostream& operator<<(ostream& out, const SpectraAlignment<U, V>& spalign){
    int iB=0;
    int iA=0;
    int nA=spalign.m_spectrumA->size();
    int nB=spalign.m_spectrumB->size();

    while(iA<nA || iB<nB){
      if(iA<nA && spalign.m_index[iA]>=0){
        for(;iB<spalign.m_index[iA]; iB++){
          cout<<"                                  \t"<<((*(spalign.spectrumB()))[iB])<<endl;
        }
        cout<<(*((*spalign.spectrumA())[iA]))<<"\t"<<(*((*spalign.spectrumB())[spalign.m_index[iA]]))<<endl;
        iB=spalign.m_index[iA]+1;
        iA++;
        continue;
      }
      if(iA>=nA){
        for(;iB<nB;iB++){
          cout<<"                                  \t"<<(*((*spalign.spectrumB())[iB]))<<endl;
        }
        break;
      }
      if(iB>=nB){
        for(;iA<nA;iA++){
          cout<<(*((*spalign.spectrumA())[iA]))<<"\t                                  "<<endl;
        }
        break;
      }
      if((*spalign.spectrumA())[iA]->moz()<(*spalign.spectrumB())[iB]->moz()){
        cout<<(*((*spalign.spectrumA())[iA]))<<"\t                                  "<<endl;
        iA++;
      }else{
        cout<<"                                  \t"<<(*((*spalign.spectrumB())[iB]))<<endl;
        iB++;
      }
    }
    return out;
  }

  }
#endif
